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Abstract 

We calculated the phase diagram of a continuous system of hard spheres loaded in a quasi- 
one dimensional bichromatic optical lattice. The wavelengths of both lattice-defining lasers were 
chosen to model an incommensurate arrangement. Densities of one particle and half a particle per 
potential well were considered. Our results can be compared directly to those of the experimental 
system [Fallani et al. PRL, 98 130404 (2007)] from which our initial parameters were taken. The 
phase diagrams for both densities are significatively different to those obtained by describing the 
same experimental setup with a Bose-Hubbard model. 
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INTRODUCTION 


An optical lattice is the result of the interference of one or several pairs of laser beams to 
produce a stationary wave, felt by a neutral atom as an external potential of the form QQ: 
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This expression, that depends on the laser wavelengths (A 3 ;,Ay,A^), dehnes a three dimen¬ 
sional lattice. 14 , 14 )K are the depths of the potential barriers felt by the atoms and can 
be modihed by changing the light intensity. In a very common arrangement, two of those 
depths are hxed to values high enough to produce a tight conhnement in two directions. We 
have then quasi-one dimensional systems well described by the continuous Hamiltonian: 
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Here, V{rij) represents the interaction between the atoms i and j, separated a distance 
Tij from each other, and m is the mass of those atoms. u:± is the harmonic frequency we 
obtain at developing the hrst two terms of the right-hand side of Eq. ([T]) around a; = 0 
and ?/ = 0, a good approximation when 14 and I 4 are deep enough. When the interaction 
between particles is modeled by a hard spheres one, a system described by this Hamiltonian 
at a density of one particle per potential well exhibits superfluid-Mott insulator transitions 
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3j for values of 14 that depend on (Uj, j^. Eq. ([2]) does not consider the system to be 
harmonically trapped in the longitudinal direction, 2 ;. That would mean to include an 
additional term of the form in the Hamiltonian 

Disorder can be introduced on this simple arrangement. Two are the most common 
experimental routes: the use of two superimposed lasers of different (and incommensurate) 


waveleng 
lattice [9 


hs on the direction 


ffl. 


or the introduction of a speckle held on the optical 
10| . In both cases, we end up having potential wells of unequal depths, with 


a quasi-periodic variation in the hrst case and a more randomly distribution of potential 
minima in the second. What the phase diagrams of those arrangements have in common 
is the presence of at least a superhuid (SF) and a Bose glass (BG) phase, depending on 
the particular experimental setup. If the average number of atoms per potential well is 
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an integer, a Mott insnlator (MI) can be also stable for some valnes of the lattice defining 
parameters [m |. 

In this work, we considered the hrst disorder option, and solved the Schrodinger eqnation 
corresponding to the continnons Hamiltonian in Eq. (j2]) with an optical lattice potential of 
the form: 


Vext{z) = siEr^ sin^ 
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where Ai and A 2 are the two wavelengths that dehne the bichromatic lattice. Here, 
and Er 2 are the recoil energies corresponding to Ai and A 2 (Er^= h^/(2mA^)). Since we 


chose to reprodnce as closely as possible the experimental setnp nsed in Ref. y, we took 
Ai = 830 nm, = 40 (in nnits of E^J, and m as the mass of a ®^Rb atom. There is a 
difference, tongh. We considered A 2 = 1067 nm instead of the 1076 nm originally nsed. This 
makes the ratio A 1 /A 2 = 7/9 ~ 0.7778 instead of the 0.77137 of the experimental data. si 
and S 2 where taken to be in the same range as the experimental data (si = 16; 0 < S 2 < 
2.5, in nnits of their respective recoil energies). To nse rational approximations to the real 
ength ratio is a common practice (see for instance Ref. Il2|, in which A 1 /A 2 = 5/7, or 


wave 


Ref. 


13l in which A 1 /A 2 = 97/126), and allowed ns to nse relatively small simnlation cells. 
The interatomic potential was described by a hard-spheres (HS) interaction: V{rij) = -|-oo 
for Tij < a and 0 otherwise, a being hxed to a single valne, that of the scattering length of 
the ^^Rb atom s |l4l| . This potential was previonsly nsed to describe bosons loaded in optical 


lattices Il5 


m- 


In most previons literatnre, the continnons Hamiltonian of Eq.(l2]) is approximated by the 
discrete Bose-Hnbbard (BH) model 20l |: 
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in which all the particles are allowed to stay only at the minima of their corresponding optical 
lattice potentials {i sites). Besides, only interactions of the type < ij > were considered, 
j being one of the two nearest neighbors (a qnasi-one dimensional system translates into 
a pnre one-dimensional Hamiltonian) to any particnlar i. 5/’(5i) is the bosonic creation 
(annihilation) operator at i, and n* represents the nnmber of nentral atoms at that same 
location. J, U and e, are parameters that can be related meaningfnlly to the emerimental 
ones ( 14 , A, and wj,) only when the optical lattice potentials are deep enongh [2, l20t]. For 
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the kind of Hamiltonians we are dealing with in this work 
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J is the hopping matrix element between sites i and j, and U models the repulsive inter¬ 
action between two atoms located at the same i site. To deduce the expression for U, it 
was supposed that the full interaction potential between atoms in the original continuous 


Hamiltonian was adequately described by a pseudopotential 


0 . 


Eqs. 


and 


were 


2l| . However, this is usually not 


deduced for monochromatic lattices, i.e., lattices dehned by a single laser of wavelength A, 
and recoil energy Er. If another laser is introduced, both expressions would have to be 
recalculated to include the effect of this second wavelength 

done 3, 3, 21-25], even though a variation in J could produce an appreciable change in 
the phase diagram. 26] 

The remaining parameter in Eq.(jl]), Aj, represents the energy offset at each potential well, 
and for a monochromatic lattice with no longitudinal conhnement is the same for each site. 
This means that could be set to zero without loss of generality. However, when disorder 
is introduced, this is no longer the case. For bichromatic lattices. A, oc cos[27r(Ai/A2)?] 
8, 3, while for random disorder, Aj follows a uniform distribution in the interval 
[-A,0] or [-A, A] pl-bsj. 

In this work, we calculate the phase diagram of the continuous Hamiltonian dehned by 
Eq. ([2]) with the optical lattice potential of Eq. ([3]). This implies that the disorder arises 
naturally at the introduction of the second laser, with no ad hoc change in the depths of the 
potential wells. The only approximation was to use a hard spheres interparticle interaction. 


a very common approach 


15 - 17 , 3 , 19 


19|. We considered two densities, one particle per 


optical lattice minimum, and half a particle per minimum, and compared our results to those 


calculated with a Bose-Hubbard model for the same systems 


2l| . However, our primary 


goal will be to compare our results to the experimental ones, when available. Since the 
lattice dehned by Ai is in most cases deeper than the dehned by A 2 {siEr.^ > S 2 ErY, to 
dehne the average density we used the number of potential minima corresponding to the 
lattice dehned by Ai. This convention will not ahect the limits of the phase diagram, since 
when SiEr.^ ~ S 2 Er.^, the system is always a superhuid, as we will see below. 
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METHOD 


The initial step to calculate the phase diagrams was to obtain the ground state of a set of 
atoms described by the Hamiltonian of Eg. (1^. T his means that we are operating a T=0, a 

131 . 21| to describe the very low temperature 


reasonable approximation commonly used |12l. 


regimes in which the experiments are performed 


100 . 


In any case, we expect our description 


should be valid for temperatures lower than the corresponding to the excitations of the 


system, which are in the range of the kHz 0. To do so, we used the diffusion Monte Carlo 
(DMC) method [^. That technique needs an initial approximation to the ground state 
wavefunction to work, the so-called trial function. That initial expression depends on the 
positions of the N atoms in the simulation cell (ri,- • • ,r 7 v). In this work, we used: 

N N N 

*I>(ri, • • ■ ,rN) = W'ilj{,Xi,yi)W(t){zi) n (7) 

i=l j=l l<m=l 

where x*, ^/j, Zi are the coordinates of the atom i and rim is the distance between atoms 
I and m. We took tfj^Xi^yi) to be a Gaussian of variance = h/{mu±), i.e., the exact 
solution to the harmonic potential that traps the particles transversally. (p{zi) was taken 
to be the numerical solution to an one-dimensional Schrodinger equation including only 
the monochromatic lattice defined by Ai. The two body correlations were modeled by: 
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a is here the scattering length of the pair of atoms. T 
constants in this equation has been detailed elsewhere 


o 

Tij > D. 


( 8 ) 


le particular way to obtain the five 

18|,ll9|. 


As indicated above, we chose the si and S 2 parameters in the same range that the ex¬ 


perimental ones of Ref. 0 The main difference is that instead of limiting ourselves to the 
case si/Er^ = 16, we considered the entire range 0 < si/Er^ < 16. For each si/Er^ we 
performed calculations in the interval 0 < S 2 /Er^ < 2.5 for a density of one atom per po¬ 
tential well and 0 < S 2 /Er^ < 2 when the average atom occupation was 0.5. In most cases, 
a simulation cell including 36 potential minima (= 18 Ai) was used. We checked that to 
increase that number to 54, 72 or even 90 did not vary the results presented in this work. 
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RESULTS 


One particle per potential well 


A Bose glass 
(gapless) phase 


is usually defined as an insulating (non-superfiuid) and finite-compressibility 
11 


ll. I23l| ■ On the other hand, a Mott insulator is a non-superfiuid gaped state. 
So, to determine whether a particular arrangement of atoms can be described as a glass, we 
should calculate the superfiuid fraction (to know if we are dealing with an insulator) and the 
compressibility of the system (to see if we have a MI or a BG). We define the compressibility, 
K, as: 

where fi is the chemical potential, 

dE E 
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Here, E is the total energy of the hard spheres arrangement loaded in the quasi-one dimen¬ 
sional optical lattice and n is the density of the system (n = iV/L, N, number of atoms; L, 
length of the simulation cell). This means that a finite k at a given n is associated with a 
continuous slope of the energy per atom versus the density at that n. 

In Fig. [U we represent E/N, in units of the recoil energy of the primary lattice versus its 
filling fraction (uAi/2), subtracting from the energy the constant contribution per atom due 
to the transversal harmonic confinement, Eho/N, in the same units. All the cases displayed 
correspond to Si/E r^ = 16 and different values for S 2 . S 2 /Er^ = 1.2 (bottom), S 2 /Er^ = 1.3 
(middle), and S 2 /Er^ = 1.4 (top). To see the influence of the size of the simulation cell in 
our results, we display there the energies corresponding to simulation cells including 36 (full 
symbols), 54 (open symbols) and 72 (crosses) potential minima. The results for bigger cells 
were similar and are not shown for simplicity. The dashed and full lines are least-squares 
fits to the simulation results for densities below and above nAi/2 = 1 when only the data 
for 36 wells are considered. The fact that all the symbols are basically on top of those fits, 
indicates that our energies are not affected by size effects. The dotted line corresponds to 
the continuation of the result in the 0.85 < nAi/2 < 1 range beyond that last density. That 
line helps us to see that for S 2 /Er^ = 1.3, the slope of the energy versus the filling fraction 
is discontinuous at nXi/2 = 1, something similar to what happens for S 2 /Er^ = 1.2. 
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FIG. 1: (Color online) Energy per particle (minus the harmonic contribution, Eho/N) for si = 
16 in units of versus the filling fraction. Symbols correspond to the simulation results for 
different values of S 2 (from bottom to top S 2 = 1.2, S 2 = 1.3 and S 2 = 1-4, all of them in units of 
Eji^). Full symbols correspond to simulation cells with 36 potential wells, open symbols to cells 
with 54 minima and crosses display the case for 72 wells. Lines correspond to least squares hts to 
the simulation data considering only the cell with 36 potential minima. 


On the other hand, when S 2 /Epi^ = 1.4 there is no such a slope change, what indicates 
that for those parameters we have a gapless phase. Using that information, and since we 
considered the smallest difference between consecutive S 2 /values to be 0.1, we estimated 
the upper limit of the Mott insulator for si/E r^ = 16 to be S 2 /Er^ = 1.3 ± 0.1. This corre¬ 
sponds to the rightmost open symbol in Fig. [2J For larger S 2 /Er^ values, we have a gapless 
arrangement that can be either a superfluid or a Bose glass. That value is in reasonable 
agreement with the experimental results of Ref. y, in which we can see a disappearance of 
the MI somewhere between S 2 /Er^ = 0.5 and S 2 /Er^ = 1.5. The remaining open symbols 
in Fig. [2] were obtained in the same way that the already described for si/Er^ = 16 and 
indicate the upper limit of stability of the Mott insulator phase for each si/Er^ value. 

We calculated the superfluid fraction (p^/ p) using a variant of the winding number tech- 
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FIG. 2: (Color online) Phase diagram for the case of an average density of one particle per potential 
well. SF, superfluid; BG, Bose glass; MI, Mott insulator. si and S 2 are the parameters defined in 
Eq.Q, in units of their respective recoil energies. See further explanation in the text. 


nique 


29|, l30|]. The formal expression employed was: 
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where Dq = fi?/2m. r indexes each Monte Carlo step of a given simulation, from the initial 
step after the Monte Carlo equilibration (r= 0) to the end of the Markov chain (in our case, 
r = 20000). D{t) is dehned as: 


0{t) = - rcM(0)f} 


( 12 ) 


Here, fcM{T) is the position of the center of mass of all the particles in the simulation cell 
at a step labeled by r. N is, as before, the number of particles. This dehnition means that 
if we represent D{t)/Dq as a function of r, for large enough r, we should get a straight line 
whose slope will give us the value of the superfluid fraction. An example of that procedure is 
displayed in Fig. [31 Each curve represents the average of 50 independent DMC simulations. 
There, S 2 /Efi^= 1.5 in all cases, while Si goes from si/Er.^ = 16 (bottom) to = 5 

(top). The solid line corresponds to the case of a perfect superfluid {D{t)/Dq = rps/p = r), 
and it is included as a reference. The error bars are displayed only for the case Si/E^^ = 
5, are fairly typical, and correspond to an standard deviation of the data produced in those 
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FIG. 3: (Color online) Diffusion of the center of mass as a function of the Monte Carlo step index, 
r. All cases correspond to S 2 /Eii^ = 1.5 and different values of si/Er^. From bottom to top: si 
= 16,12,10,8,5. The solid line gives the slope corresponding to Ps/p = 1- The error bars are only 
displayed in one case for simplicity. 


50 independent simnlations. What we can see is that the slope of D[t)/Dq for large valnes 
of r (> 10000) increases from zero = 16,12), to valnes that imply the existence of 

snperflnids [si/Er^ = 10,8,5). 

From the slopes of cnrves similar to those in Fig. [3l and with the help of Eq. ffTTj) . we 
can dednce the snperflnid fraction for the set of parameters we are interested in. Fig. 0] 
shows those results for different si/Er^'’s as a function of S 2 /Eji^. The errors bars for each 
Ps/p were calculated by the same method as in Ref. la least-squares fitting procedures were 
performed for the average values of D{t)/Dq as a function of r, and for the same values after 
the subtraction (or addition) of the corresponding limits given by the error bars of that ratio. 
This allows us to obtain a range of most probable slopes, and from that, a superfluid fraction 
interval. What we observe is that, except for the si/Er^ = 5 case, there is a range of S 2 /Er^ 
values above which the system is an insulator, and below which we have a superfluid. We 
consider an arrangement to be non-superfluid when the error bar for the superfluid fraction 
includes the value ps/p = 0. So, for instance, when si/Er^ = 12 (full squares in Fig. Hj), 
the last insulator corresponds to S 2 /Er^= 1.6, since in that case Ps/P= 0.01 ± 0.01. This 
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FIG. 4: (Color online) Superfluid fraction as a function of S 2 (in units of Er^) for different si values 
(in units of Er^). Open squares, si/Er^ = 16; full squares, si/Er^ = 12; open circles, si/Er^ = 
8; full circles, si/Er^= 5. The error bars are only displayed in the last two cases for simplicity, 
the ones not shown being of the same size {A{ps/p) 0.01). 

translates into the point Si/Er^= 12, S 2 /ER 2 = 1.6 in Fig. |2l For larger S 2 /Er^^s, we have 
a snperflnid, and when 1.3 < S 2 /Er^ < 1.6, the phase is a gapless insnlator. i.e. a Bose 
glass. The remaining of the fnll symbols in Fig. [2] have been obtained in a similar way for 
other si/Er^ valnes. The error bars displayed there have the same origin as the ones for 
the MI-BG limit. The fnll symbols in Fig. [2] are related to the snperflnid behavior of the 
system (see below). 

Half a particle per potential well 

Exactly the same calcnlations can be performed for an average density of half a particle 
per potential well. The only difference is that in the phase diagram we will have only 
snperflnid and Bose glass phases, since a MI is only possible for a hlling fraction of one. 
Fig. [5] is the connterpart of Fig. IHfor this lower density. There, we displayed the snperflnid 
fraction as a fnnction of S 2 /Er^ for si/E r^ = 8 (npper curve) and si/Er^ = 10 (lower curve), 
both sets of calculations including their corresponding error bars for each point. We can see 
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that for si /= 8, the system is always a superfluid, irrespectively of the value of S 2 /Er^. 
On the other hand, for si/Er^ = 10, an insulator (Bose glass) appears for S 2 /Er^ > 0.9. 
This last value is represented in Fig. Elas the leftmost full symbol. Additional calculations 
indicate that for this density no Bose glass appears when si/Er^ < 8, at least for S 2 /Er^^s 
in the range displayed in Fig. O For larger si/Er^^s, the limits between the superfluid and 
insulator phases are also shown as full symbols. The upper dashed line is a least-squares £t 
of those simulation results to a function of the form S 2 (si) = a/(si — 8)^ (a, b, constants, S 2 
and Si in units of their respective recoil energies), i.e., to a function with an asymptote for 
Si/Er-^ = 8 (vertical dotted line), and it is intended simply as a guide-to-the-eye. 

On the other hand, the open symbols in Fig. O correspond to the values obtained for a 


Bose-Hubbard Hamiltonian that models the same system 


2l| . For instance, the rig 


it most 
{U/J 


open symbol of Fig. [6] correspond to the largest U/J value considered in Ref. [21 
= 10), for which the limit between the superfluid and the Bose glass phases is located at 
S 2 Er^/J = 4. Using Eqs. [5] and [6] above, we have that this translates into si/Er^ = 6.1 and 
S 2 /ER 2 = 0.4. A similar procedure was used to obtain the remaining open symbols shown 
there, the only cases whose transformations give coordinates in the S 2 /Er^ range displayed 
in Fig. [6l The dashed line on top of them is again a least-squares £t to be used only as 
a guide. What we observe is that both in our simulations and those of Ref. l2l| predict a 
Bose glass for high enough S 2 /Er^ and Si/E r^, even though the Bose-Hubbard Hamiltonian 
underestimates the values of those parameters needed to create an insulator. 


DISCUSSION 

In this work, we described the behavior of a set of neutral atoms loaded in a quasi- 
one dimensional bichromatic optical lattice. To do so, we solved the Schrodinger equation 
corresponding to the continuous Hamiltonian in Eq. ([2]), with an external potential dehned 
by Eq. [3l The only approximation employed, besides the use of a stochastic method to 
obtain the ground state, was the description of the two-body interaction by a hard-spheres 
potential. Within this framework, we obtained the phase diagrams of Figs. |2] and [6l as 
functions of the real parameters at two different atom densities. This makes straightforward 
the comparison between theory and experiment. In particular, when si/E r^ = 16, our value 
for the MI-BG limit was S2/Er^=1.?) ± 0.1, a value within the experimental range (0.5 
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FIG. 5: (Color online) Same as Fig. HJ but for a density of half a particle per potential well. Full 
squares, si/En^ = 10; open squares, si/Er^ = 8. 
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FIG. 6: (Color online) Same as in Fig. O but for an average density of half a particle per potential 
well. The full symbols correspond to our simulation results while the open symbols are from Ref. 
2ll . All the lines displayed are guides-to-the-eye. See further explanation in the text. 
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< S 2 /-ER 2 <1.5) j^. 

However, the standard approximation to describe an optical lattice is the nse of a Bose- 
Hnbbard model. In a classic case of pntting the cart before the horse, the fact that a set of 
atoms loaded in an monochromatic optical lattice is a very good experimental realization of 
a BH Hamiltonian for some valnes of the lattice defining parameters js, has been nsed 
as an excnse to describe any optical lattice in any circnmstances by this model. The main 
advantage of this approach is that the BH model is relatively easy to nse and its behavior 
is well nnderstood; the main problem being that it is only an approximation that breaks 
down even for monochromatic lattices with shallow potential wells, as it has been established 
both experimentally M| and theoretically jd, 18| for the description of the snperflnid-Mott 
insnlator transition. 

In addition, the stndy of a dichromatic lattice implies corrections for both the J and U 
parameters that wonld inclnde the effect of the second standing wave and that there are 
nsnally not considered. This means that even in those cases in which the primary lattice 
is deep enongh to be well described by a BH Hamiltonian, the description of the lattice 
is not exact. That componnded error wonld be the canse of the differences between the 
phase diagram for half a particle per potential well calcnlated with a continons and a BH 
Hamiltonian shown in Fig. |6l Unfortnnately, to onr knowledge, there are no experimental 
data to compare those diagrams to. We can only say that the stability zone of the snperflnid 
phase is rednced with respect to the case of a continnons Hamiltonian, the same general trend 
fonnd for monochromatic lattices j4|. 

The sitnation is similar for a density of one particle per potential well. As indicated 
above, for s\/Er^ = 16, there is a reasonable agreement between the position of the MI- 
BG transition predicted from onr calcnlations and the experiment. Regretably, we cannot 
compare any BH predictions to that valne, since Si = 16 translates (nsing Eqs. ([5]) and ([6]), 
with the scattering length taken from Ref. Q) into U/J 130, a larger valne than the 
ones considered in the 


the other hand, the npper stability limit 


Lguii udJivcii 

literatnre [12, 21,l22[. On 


for the MI phase in Ref. l2l| for U/J = 10 (larger U/J nsed) corresponds to S 2 Er^/J ~ 10. 
This means Si/Er^ ~ 6.1 and S 2 /E/J 2 ~ 1-1, to be compared to Si/Eji^ = 6, S 2 /EJI 2 = 0.6 
± 0.2, from Fig. |6l Thns, we have that a BH description prodnces a shrinking of the SF 
stability zone with respect to the case of a continnons model, as in the case of the lower 
density described above. This is corroborated by the onset of the SF-MI transition for S 2 
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= 0. The value in this work {si/Er^ = 5 ± 1) is larger than the one deduced in Ref. [21 
{si/Er^ ~ 3). 

To conclude, we should make a last remark about the general form of the phase diagram 
displayed in Fig. [2J There, we can see that for some values of Si/Er^ [si/Er^ > 10), 
an increasing in S 2 /Er^ makes the system change from a MI to a BG and further up, to 
a superfluid. However, when si/Er^ is smaller, there is only a single transition, MI 
SF. This is in apparent contradiction to some previous theoretical results [n|, Q for a 
Bose-Hubbard model, that indicate that between a MI and a SF phases, there should be 
always a BG ’’sandwich”. At least part of the discrepancy can be explained by the fact 
that those works correspond to systems with random site disorder, which behave differently 
than arrangements with quasiperiodic dichromatic lattices. In any case, the results of Fig. 
|2] indicate that if that intermediate BG exists for si/Er^ < 10, the width of its stability 
interval should be very narrow. 
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